Harmonic inversion helps to beat time-energy uncertainty relations. 



Zbyszek P. Karkuszewski 
Institute of Physics, Jagiellonian University, ul. Reymonta 4, 30-059 Cracow, Poland 

(Dated: February 1, 2008) 

It is impossible to obtain accurate frequencies from time signals of a very short duration. This 
is a common believe among contemporary physicists. Here I present a practical way of extracting 
energies to a high precision from very short time signals produced by a quantum system. The product 
of time span of the signal and the precision of found energies is well bellow the limit imposed by 
the time-energy uncertainty relation. 
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INTRODUCTION 

Many methods of obtaining frequencies or energies 
from time signals have limited resolution. This limita- 
tion is often expressed in the form of a time-energy un- 
certainty relation. For instance, the width AE of spectral 
lines in atomic spectroscopy is determined by life time r 
of the excited atoms 

AEt = h. 

Another example is a Fourier transform applied to a dis- 
crete time signal in the interval of length T. The grid 
step in a frequency domain and thus the resolution of 
the method is given by Aw = 2n/T, which in turn yields 

AujT = 2tt. 

All relations of this kind may give a wrong impression 
that something more fundamental is the source of the 
uncertainties - a principle that one cannot measure en- 
ergies to an arbitrary precision in a very short time. Such 
an intuition is common if the system subjected to energy 
measurement is quantum. 

Uncertainty principles in quantum mechanics are a 
mathematical consequence of the following theorem. 
Two self-adjoint operators A and B defined on the same 
Hilbert space necessarily obey the relation 



AAAB > -I 
~ 2 1 



\A.B] 



(1) 



where [...,...] is a commutator, (...) stands for quantum 
average in a given state from the domain of A and B, 
and (AZ) 2 = (Z 2 ) — (Z) 2 . In the case of position x and 
momentum p operators gives 



AxAp > 



(2) 



Notice that Ax and Ap are completely determined by 
the state of the system being measured and have nothing 
to do with an accuracy of the measuring apparatus. lf2"| 
holds even if this accuracy is infinite. The correct ap- 
proach to the quantum time-energy uncertainties can be 
found in Q. The theorem does not apply to the case 
of time and energy, simply because time is not an oper- 
ator, it is just a parameter in quantum mechanics. One 



gains nothing forcing the idea that a parameter is a spe- 
cial kind operator because JQ) gives zeros on both sides 
and does not provide grounds for existence of any time- 
energy uncertainty relation. Indeed, it has been shown 
how to precisely measure energy in an arbitrarily short 
time . Some methods of the energy measurement have 
their own limitations |2j reflected in loss of accuracy when 
applied for very short time signals. Here I use the method 
free of such inconveniences. 

Let us setup some general issues before getting to the 
details of the method. Suppose that a continuous signal 
c(t) is given in a finite time interval 



c(t) 
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for < t < T, 



(3) 



where K is finite, dk is a real amplitude and cuk is a 
real frequency. Complex frequencies will be considered 
in future. The task is to find unknown amplitudes and 
frequencies of c(t). What one can see is that c(t) is an 
analytic function of time t and, as such, can be uniquely 
extended beyond the interval (0, T). This means that, in 
principle, even for tiny T it is possible to get all ampli- 
tudes and frequencies to a high precision from J2J. Unfor- 
tunately it seems difficult to solve this nonlinear problem 
analytically and numerical methods cannot handle con- 
tinuous signals due to infinite number of data points. One 
way of getting around this problem is to take only finite 
number of points at cost of loss of uniqueness of the ex- 
tension. However, as will be shown later, this drawback 
is usually not severe for practical purposes. From now on 
I shall assume that the signal c(i) is known only at N + 1 
equidistant time points t n — nSt for n — 0,...,N and 
tjsi = T. 10 can be rewritten as a set of A + 1 equations 
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(4) 



where c n = c(n8t). This set has 2K real unknowns (K 
amplitudes and K frequencies) so the number of equa- 
tions N + 1 has to be equal or greater than K. This 
condition would have opened the possibility of existence 
of a unique solution if the equations were linear in lu^ 
and dk- It is not the case here. There will always be infi- 
nite number of solutions to (@J if the set is self-consistent 
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and no solutions otherwise. For reasons explained in the 
next sections it will be required that N > K. Notice 
that discrete Fourier transform method assumes the grid 
of K equidistant frequencies and solves Q only for d k as 
a linear set of equations. The more challenging task of 
solving 10} for amplitudes and frequencies is performed 
by, so called, harmonic inversion method. 



INTERVALS OF UNIQUE SOLUTIONS 



Namely |$ n ) = J^feLi a fc u fcl u fc), where a k stands for a 
time independent complex number. JfjJ) in the new form 
appears as 



The two functions e M and e tuJ have the same 
values at equidistant time points nSt, n = 0, N if 

Ll> = Ul + 

not 

The integer number I has to be chosen in such a way 
that the fraction l/n is integer as well. This implies that 
I = mN\, where m is integer and 



2ttN\ 



(5) 



The time discretization of the signal c(t) , the transition 
from to Q , happens at cost of the loss of uniqueness 
of solutions. However, as can be seen from solutions 
are unique in finite intervals in frequency domain. Fre- 
quency u> is unique in the interval 

27riViV! 2nNN\\ 



T 



T 



which can be made large by increasing N and/or decreas- 
ing T. Conversely, if one posses prior knowledge about 
the range of frequencies in a signal it is straightforward 
to design a sampling step St to extract all frequencies in 
a unique way. 



HARMONIC INVERSION 

The key idea behind the harmonic inversion method 
is to replace the original nonlinear problem Q with an 
eigenvalue problem of an operator, as in J^J. The pre- 
sentation of the method in this section follows that in 

nt 

Lets start with a normalized quantum state |$o)- The 
evolution of this state is generated by unitary evolution 
operator U(8t) and |$„) = U n (5t)\<f> ) . The states |$„) 
are so important that they were given a name of Krylov 
states. For every time signal in there exists such an 
evolution by time St operator U that the signal can be 
viewed as an autocorrelation function 



($ |$ r 



(6) 



Krylov states in ijfJJl can be written in the orthonormal 
basis of eigenvectors of U defined by U\uk) = u k \uk)- 
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Comparing JJJ) and one can see that eigenvalues 
Uk = e~ lulkSt and |afc| 2 = dk- It is enough to find eigen- 
values of U in order to obtain all frequencies u>k in the sig- 
nal c„. These frequencies multiplied by Ti can be viewed 
as eigenenergies of a hypothetical Hamiltonian governing 
the evolution of a system in the initial state |$o)- Of 
course, one can start from the autocorrelation function 
© and find the eigenenergies of the real system. 

It turns out that the matrix elements of U in the basis 
of Krylov states can be expressed in terms of c„ alone, 
i.e. without explicit reference to the states 



-»+l 



$n) = C 



(8) 



where i,j — 0, ...,N — 1. The negative indices of c 
in the equation above introduce no complication since 
c_„ = c*. To obtain K eigenvalues the dimension N of 
the matrix U must be equal or greater than K. It means 
that the number of complex signal points c„ required by 
the method exceeds the half of the number of unknowns. 
The Krylov vectors are normalized but not orthog- 
onal. Thus the eigenequation for matrix U in the Krylov 
representation © takes the form 



U\uk) = u k S\u k ) 



(9) 



where S is a matrix of scalar products Sij = ($i|<&j) = 
Cj-i with i, j — 0, N — 1. 

The harmonic inversion method consists of two stages. 
First is to solve generalized eigenvalue problem @ , where 
all matrix elements are expressed in terms of c n , in order 
to get all frequencies. Second, when all frequencies in (@J 
are known, it is enough to solve linear set of equations for 
the amplitudes d k ■ There are some practical difficulties in 
proceeding with the first stage. For instance, numerical 
algorithms fail in finding the eigenvalues u k if the signal 
duration T is small and the number of frequencies K > 4. 
The second stage is straightforward and will be skipped 
in this work. It will be assumed that all real amplitudes 
dk are equal and the signal is normalized to unity, i.e. 
c = 1. 

In the next section we will see that the crucial role in 
the numerical approach to JHJ is played by the smallest 
positive eigenvalue of S. 



PROPERTIES OF MATRIX S 

One way of dealing with the generalized eigenvalue 
problem is to reduce it to the ordinary eigenvalue prob- 
lem by multiplying both sides of © by the inverse of 
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the Hermitian Toeplitz matrix S. However, matrix S is 
singular unless N — K and the inverse does not exist. 
Indeed, recall that every vector |$ n ) can be decomposed 
into a superposition of K linearly independent eigenvec- 
tors of U. Therefore the rank of the matrix S is K. 
Additionally, those nonzero eigenvalues are positive be- 
cause a scalar product (y>\tp) > for any nonzero state 
\ip) in a Hubert space. Suppose that columns of a N x K 
matrix G are represented by eigenvectors of matrix S 
to the positive eigenvalues. In order to restrict © to 
the range of S we define W = G^UG, S' = G^SG and 
\u' k ) = G~ x \u k ). S' is diagonal and positive defined thus 
the ordinary eigenvalue problem converts to 



S'- 1 U'\u' k )=u k \u' k ). 



(10) 



Construction of matrix G is possible only if one can ex- 
tract all positive eigenvalues of S. This task becomes 
hopeless if the smallest eigenvalue of S cannot be dis- 
tinguished from zero due to limited accuracy. The ex- 
act analytical expression for the smallest eigenvalue of a 
Toeplitz matrices is not known yet. Here is the approx- 
imate formula for the magnitude of the smallest eigen- 
value for the short time span T of the signal 



A, 
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(11) 



f2 is the magnitude of the greatest (to the absolute value) 
frequency in the signal. For the derivation of i|ll[l see 
the directions in the Appendix. The formula is valid 
only if the shortness condition is fulfilled, TVl <C 1 . The 
expression in square brackets in 1)11(1 is smaller than 1 and 
the Xmin decreases exponentially fast with increasing K . 
For the estimated magnitude of the function a(w\, luk) 
see Fig^ 

The eigenvalue X m in is additionally corrupted by inac- 
curacies of the signal. The accuracy analysis is presented 
in the next section. 



IMPACT OF NOISE 

From now on it will be assumed that the accurate sig- 
nal c„ is perturbed by noise r\ n and the new signal 

Cn C n -\- Tj n . 

The noise is limited, r\ n G (—r) max ,r] max ) for all n and 
Tjmax > 0. Matrix S, which is constructed using c n in- 
stead of c n (see the text under ©), remains Hermitian. 
It can be shown || that the smallest eigenvalue of S must 
obey the inequality 

|A m j n — A m .;„| < 2Nri max + 0(rj max ). 

If X min is to be distinguished from perturbed zero eigen- 
values the following lower bound on X m i n arises 
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FIG. 1: Statistical properties of a(uii, ujf). For every K 
100 sequences of K frequencies were randomly generated. 
Each frequency has been uniformly drawn from the inter- 
val (0.5, 1.0). For each sequence a(ui, ...,Uk) was computed. 
Circles denote the most probable value of a(oj\, ...,ojk) and 
vertical bars embrace 90% of sequences. 



This diagnostic condition combined with ((11(1 provides 
the limits of applicability of the harmonic inversion 
method for short signals. 

Tracking the inaccuracy propagation when solving 1(10(1 
one arrives at the surprising at first sight "certainty re- 
lation" on perturbed frequencies u>k 



\oj k -0Jk\T < 



2KN 2 



X r , 



(13) 



Xmin > 4A?7„ 



(12) 



The error concerning frequencies of short signals can be 
made arbitrarily small by reducing the noise amplitude! 
Higher order terms in r\ max were dropped in ((13(1 . For 
time signals of short duration the presence of X m in in the 
denominator in RHS of ((13(1 is unfavorable. As stated in 
lfTT|l X m in rapidly goes to zero as K increases. Therefore, 
to extract many accurate frequencies from a short signal 
the very low level of the noise r\ m ax will be required. 

In numerical experiments the noise is caused by finite 
precision of floating point numbers. As an example, the 
harmonic inversion method was used on a signal with 
K = 10 frequencies drawn from an interval (0.5,1.0), 
sampled at N = 14 points with T = 0.01 and using 85 
digits precision i.e. i] max — 10~ 84 . The results are shown 
in Table Notice that the Fourier algorithm applied to 
this signal would give the grid step in frequency domain 
200-7T which is several orders of magnitude greater than 
the error in Tabled In the example above X m i n — 4.07 x 

io- 78 . 



DISADVANTAGES OF THE METHOD 

Harmonic inversion method does pretty well when ap- 
plied to short time signals provided the level of inaccu- 
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TABLE I: Numerical example. All amplitudes d k were equal. 







\u>k - 


0.50415486481506 


0.50415486481507 


0.00000000000001 


0.51315149664879 


0.51315149664880 


0.00000000000002 


0.66526505816728 


0.66526505819813 


0.00000000003085 


0.71068158128764 


0.71068158894354 


0.00000000765589 


0.73253390251193 


73253391404702 


00000001153508 


0.75819833122659 


0.75819832230088 


0.00000000892572 


0.79694000270683 


0.79694000183578 


0.00000000087105 


0.85043252643663 


0.85043252642270 


0.00000000001394 


0.88220469909720 


0.88220469909823 


0.00000000000103 


0.93358750757761 


0.93358750757763 


0.00000000000001 



racies in the input data is very small. 

In principle this level can be kept very low in mea- 
surement of quantum systems. For instance, the mea- 
surement of quantum autocorrelation function squared 
(probability of finding a system in the initial state at 
different times of its evolution) can be arbitrarily pre- 
cise if done simultaneously on many copies of the system. 
The inaccuracy decreases as an inverse square root of the 
number of the copies. This function involves (K— l)K/2 
frequencies rather than if as in JJjJ. 

In experimental practice, however, inaccuracies are big 
and the method might be useful only for signals with 
small number of frequencies. 

Even if desired accuracy is provided another problem 
may show up. The harmonic inversion method as de- 
scribed here involves diagonalization of N x N matrix S 
and K xK matrix S^ 1 !! . S matrix is Hcrmitian Toeplitz, 
the matrix elements are constant along diagonals, and it 
is enough to store only one row of S. Solving a linear set 
of equations with Toeplitz matrix requires oc N log2 (N) 
operations using superfast algorithms 0. Perhaps di- 
agonalization of Toeplitz matrices can be also speed up 
bellow oc A^ 3 operations. The complexity of the second 
diagonalization is oc K 3 operations but here no eigen- 
vectors need to be computed. Remember that all these 
operations has to be done with high precision and are, 
therefore, relatively slow. 

The diagonalization related difficulties has been over- 
come for long signals 0, @- F° r such problems the 
frequency domain was effectively divided into smaller in- 
tervals containing fewer frequencies and the harmonic in- 
version, known there as filter diagonalization, was applied 
to one interval at the time. Splitting the frequency do- 
main introduces however its own inaccuracies into the 
computed frequencies of the signal and filter diagonaliza- 
tion methods are said to obey their own time-frequency 
uncertainty relations. The total duration of the signal is 
limited from bellow by local density of frequencies 0] or 
by minimal and average frequency distance 0. It is not 
clear yet whether these restrictions can be invalidated by 



improved precision of calculations. 

SUMMARY 

In this work I presented the application of the har- 
monic inversion method to short time signals. I have 
demonstrated that the method has no fundamental lim- 
itations regarding the length of the signals. It works for 
arbitrarily short signals if sufficient accuracy of the in- 
put data is provided. In particular, it has been shown 
that it is possible to extract energies from the short au- 
tocorrelation function generated by a quantum system. 
The method can be also applied to short autocorrelation 
functions squared which have a clear experimental inter- 
pretation. Its measurement can be, in principle, carried 
out in arbitrarily short interval of time and still all ener- 
gies involved in its evolution can be extracted to a desired 
accuracy. 
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APPENDIX: ESTIMATE OF \ min 

The direct calculation of A m , n is hard. What can 
be done without much effort is the prediction of the 
general form of the expression for this quantity. The 
first step is the observation that since each matrix el- 
ement of S is a sum of K exponents as in Q then 
S can be written as a sum of K matrices S = S\ + 
52 + ... + Sk, where each matrix Sk depends only on 
one frequency uj k . Matrix S k is Hermitian and has only 
one nonzero eigenvalue N and corresponding eigenvector 
\s k ) = [1, e^"" 5 * , e luJk2St , ...,e iuJk( - N -^ St ]/y/N. Eigenvec- 
tors \sk) form a nonorthogonal basis. Lets construct a 
K x N matrix R by filing its rows with vectors \sk)- 
Eigenvalues of RR) multiplied by N give the all nonzero 
eigenvalues of S. The characteristic equation for eigen- 
values of RR) 

a K {\/N) K + ... + ai(X/N) + dct(RR f ) = 

can be expanded at A = to the first order to give an 
estimate for the smallest eigenvalue 

Xmin ^ det(itRt/jjQ 
NK ~ oi 

The matrix RR) above is conveniently divided by K to 
normalize its trace to unity. For N — K matrix R be- 
comes square and det(i?) oc (VlT) k( - k ^/ 2 for fiT < 1. 
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Therefore, det^RE)/ K) w (ciftT)^- 1 ). The param- 
eter a\ is given by sum of minors of RR) / K and a\ f=a 
(c 2 £!T) < '' ft: ~ 1 - ) ( /f ~ 2 ' ) . Ci and c 2 are functions of normalized 
frequencies Wfe/O. Similar approximations are obtained if 
N > K. All calculations are summarized by the estimate 
CD- 
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